{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第七章　随机服务系统\n",
    "\n",
    "\n",
    "作者：[王何宇](http://person.zju.edu.cn/wangheyu)\n",
    "\n",
    "[浙江大学数学科学学院](http://www.math.zju.edu.cn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "排队是现代社会的一种常见现象。当一种服务系统的接待速度赶不上客户的到达频率时，就会出现排队。排队不光在人类社会出现，在很多工业设计，如 CPU 处理命令序列，动物大规模迁移等，也会出现排队。为直观期间，我们下面默认在一个针对顾客的服务系统中讨论排队模型。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 基本问题\n",
    "在一个排队系统中，由于\n",
    "\n",
    "1. 顾客到达时间和数量的不确定；\n",
    "2. 每一个顾客接受服务的时间不确定。\n",
    "\n",
    "所以对顾客而言，“我还要等多久”；对服务系统而言，“今天应该开几桌服务”，都是永恒的焦虑。这也构成了随机服务系统研究的主要内容：\n",
    "\n",
    "1. 系统的特征分布，如顾客到达、排队时间、队伍长度等随机变量的分布；\n",
    "2. 系统优化，如何配置服务参数和资源，使得在**有限**条件下达到某种最优，比如最高服务效率，或者最大服务量，等等。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 随机服务系统的组成\n",
    "\n",
    "一个随机服务系统，一般分为三部分考虑：\n",
    "\n",
    "1. 输入过程，即顾客达到的特征，主要有：\n",
    "    1. 顾客源是否有限；\n",
    "    2. 同时到达的顾客数是单个还是成批（是否并发）；\n",
    "    3. 顾客是以确定的频率到达，还是随机分布；\n",
    "    4. 顾客间的到达关系是否独立？（一个顾客是否会引发更多的顾客）\n",
    "    5. 平稳性，顾客到达时间分布是否与时间无关？\n",
    "2. 排队规则：\n",
    "    1. 损失制度，若顾客不能立即得到服务就会离去。典型的案例如停车场；\n",
    "    2. 等待制度，不能立即得到服务的顾客开始排队，具体又分为：\n",
    "        1. 先到先服务（FCFS）；\n",
    "        2. 后到先服务（LCFS）；（考虑模型问题...）\n",
    "        3. 随机服务；\n",
    "        4. 优先权服务；（VIP客户？）\n",
    "    3. 混合制，即有限队长，当顾客队列长度超过限制时，不再接受新的排队。\n",
    "3. 服务机构设置：\n",
    "    1. 服务台是否有多个并行；\n",
    "    2. 服务台一次接受单个顾客还是一批顾客；\n",
    "    3. 服务时间是随机的还是确定的；\n",
    "    4. 若服务时间服从分布，该分布是否平稳。（不随时间变化）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 服务系统模型的符号约定\n",
    "\n",
    "该符号约定由 X/Y/Z/A/B/C 组成，其中：\n",
    "\n",
    "+ X - 顾客到达时间间隔的分布，主要有：\n",
    "    * M - 指数分布；\n",
    "    * D - 确定情形；\n",
    "    * Ek - $k$ 阶埃尔朗分布；\n",
    "    * GI - 互相独立的随机分布；\n",
    "    * G - 一般随机分布；\n",
    "+ Y - 服务时间分布，也用上述记号表达不同的分布；\n",
    "+ Z - 并行服务台数量；\n",
    "+ A - 排队系统的最大容量，默认为 $\\infty$；\n",
    "+ B - 顾客源数量，默认为 $\\infty$；\n",
    "+ C - 排队规则，默认为 FCFS。\n",
    "\n",
    "所以，M/M/1 就表示 M/M/1/$\\infty$/$\\infty$/FCFS，即顾客到达和服务时间均服从指数分布，只有一个服务台，队列长度和客源无限，先到先服务，这大致是最常见的情形。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在一个随机服务系统中，我们关心的是系统中正在被服务的顾客数量 $n$，以及它在 $t$ 时刻的概率 $P_n(t)$，即 $t$ 时刻，系统中有 $n$ 个顾客的概率，也即暂态概率。这可以通过求解差分方程（离散时间）或着微分方程（连续时间）得到。自然地，对应稳态问题：\n",
    "$$\n",
    "\\lim_{t \\to \\infty} P_n(t) = p_n,\n",
    "$$\n",
    "即求稳态概率分布（顾客数量分布达到平衡，和时间无关）也是研究的对象。而评估一个系统好坏，一般采用如下指标：\n",
    "\n",
    "1. 平均顾客数 $L_s$；\n",
    "2. 平均队列长度 $L_q$；\n",
    "3. 平均逗留时间 $W_s$；(含被服务时间)\n",
    "4. 平均等待时间 $W_q$；\n",
    "5. 平均忙期 $T_b$，即服务系统两个空闲状态出现的时间间隔。\n",
    "\n",
    "计算机模拟提供了估算以上期望指标的一种方法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 基本概念\n",
    "\n",
    "**柏松过程** 设 $N(t)$ 表示在时间 $[0, t]$ 内到达的顾客数，以 $P_n(t_1, t_2)$ 表示在 $[t_1, t_2)$ 内有 $n$ 个顾客到达的概率，则顾客到达过程服从柏松分布指：\n",
    "$$\n",
    "P_n(t_1, t_2) = P\\{N(t_2) - N(t_1) = n\\} = \\frac{\\lambda(t_2 - t_1)^n}{n!} e^{-\\lambda(t_2 - t_1)}.\n",
    "$$\n",
    "当 $\\Delta t = t_2 - t_1 \\to 0$ 时，有\n",
    "$$\n",
    "P_1(t, t + \\Delta t) = \\lambda \\Delta t + o(\\Delta t)\n",
    "$$\n",
    "这里可以理解为做了欧拉-麦克劳林级数展开，所以可以认为是一个与 $t$ 无关的线性关系，与 $\\Delta t$ 成正比。或者说，从 $t$ 时刻，到 $t + \\Delta t$ 时刻出现一个新顾客的概率是稳定的常数，与历史信息无关。若将时间看作是离散的\n",
    "$$\n",
    "t_1 < t_2 < \\cdots < t_k,\n",
    "$$\n",
    "则这一特征符合马尔可夫链定义：\n",
    "$$\n",
    "P\\{N(t_k) = n \\left| N(t_1) = n_1, N(t_2) = n_2, \\cdots, N(t_{k - 1}) = n_{k - 1}\\right.\\} = P\\{N(t_k) = n \\left| n_{k - 1}\\right.\\},\n",
    "$$\n",
    "只和上一个时刻有关。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "另一方面，\n",
    "$$\n",
    "\\sum_{n = 2}^\\infty P_n(t, t + \\Delta t) = \\sum_{n = 2}^\\infty \\frac{\\lambda (\\Delta t)^n}{n!}e^{-\\lambda \\Delta t} = o(\\Delta t)\n",
    "$$\n",
    "注意这里 $\\Delta t$ 幂在 $2$ 或以上。顾对 $\\Delta t \\to 0$，可以忽略 $n = 2, 3, \\cdots$ 的情况。（对充分小的间隔，只要考虑来一个顾客的情况。）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "综合考虑，有\n",
    "$$\n",
    "P_0 + P_1 + P_{n \\geq 2} = 1,\n",
    "$$\n",
    "故在 $[t, t + \\Delta t)$ 内没有顾客到达的概率为\n",
    "$$\n",
    "P_0(t, t + \\Delta t) = 1 - \\lambda \\Delta t + o(\\Delta t) \\approx 1 - \\lambda \\Delta t.\n",
    "$$\n",
    "此外，令 $t_1 = 0$，$t_2 = t$，则有\n",
    "$$\n",
    "P_n(t) = P_n(0, t).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**指数分布** 设顾客到达时间间隔的概率分布是 $F_T(t)$，其中 $T$ 是两批顾客相继到达时间的间隔，即\n",
    "$$\n",
    "F_T(t) = P\\{T \\leq t\\}.\n",
    "$$\n",
    "换言之，等价于 $\\Delta t = t$ 出现了一个及以上的顾客的概率。再由之前泊松分布性质，$\\Delta t = t$ 等价于在 $[0, t)$ 区间至少有一个顾客到达：\n",
    "$$\n",
    "F_T(t) = 1 - P_0(t) = 1 - e^{-\\lambda t}, t > 0.\n",
    "$$\n",
    "其中\n",
    "$$\n",
    "P_0(t) = e^{-\\lambda t}\n",
    "$$\n",
    "是这段时间内一个顾客也没有的概率。由此\n",
    "$$\n",
    "f_T(t) = \\frac{d F_T}{d t} = \\lambda e^{-\\lambda t}, t > 0,\n",
    "$$\n",
    "也即 $T$ 服从指数分布，因此\n",
    "$$\n",
    "\\mathrm{E}[T] = \\frac{1}{\\lambda}, \\mathrm{var} (T) = \\frac{1}{\\lambda^2}.\n",
    "$$\n",
    "期望给出了 $\\lambda$ 的实际意义：$\\lambda$ 是在单位时间到达顾客的平均数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**平均服务时间** 一般也假设它服从指数分布，即\n",
    "$$\n",
    "f_{H}(t) = \\mu e^{\\mu t},\n",
    "$$\n",
    "其中 $\\mu$ 表示单位时间内被服务的顾客数，即平均服务率。$H$ 表示一个持续服务中的系统中相继被服务的两顾客**离开**时间的间隔。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "显然，一个服务系统的守恒性（服务中的顾客数量保持不变）本质上由 $\\lambda$ 和 $\\mu$ 决定，其中\n",
    "$$\n",
    "\\rho = \\frac{\\lambda}{\\mu}\n",
    "$$\n",
    "称为服务强度，若 $\\rho \\geq 1$，则队列会无限增加，永远达不到稳态（服务资源不够）。因此，我们一般要求 $\\rho < 1$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**M/M/1 模型** 最典型的一种服务模型，在数学上也是可解的。考虑在 $[t, t + \\Delta t)$ 内顾客数 $n$ 的变化。首先，我们忽略全部高阶项，也即 $\\Delta t$ 时刻内到达/离开两个或两个以上顾客的情况，那么剩下所需要考虑的就只有：\n",
    "\n",
    "$$\n",
    "\\begin{array}{llllll}\n",
    "\\mbox{编号}&\\mbox{原顾客数}&\\mbox{到达}&\\mbox{离开}&\\mbox{现顾客数}&\\mbox{发生概率}\\\\\n",
    "(1)&n&0&0&n&(1 - \\lambda \\Delta t)(1 - \\mu \\Delta t)\\\\\n",
    "(2)&n&1&0&n + 1&\\lambda \\Delta t(1 - \\mu \\Delta t)\\\\\n",
    "(3)&n&0&1&n - 1&(1 - \\lambda \\Delta t)\\mu \\Delta t\\\\\n",
    "(4)&n&1&1&n&\\lambda \\Delta t \\cdot \\mu \\Delta t\n",
    "\\end{array}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意这四种情况不会同时发生，故\n",
    "$$\n",
    "\\begin{array}{rcl}\n",
    "P_n(t + \\Delta t) &=& P_n(t)(1 - \\lambda \\Delta t)(1 - \\mu \\Delta t) \\\\\n",
    "&+& P_{n - 1}(t) \\lambda \\Delta t (1 - \\mu \\Delta t) \\\\\n",
    "&+& P_{n + 1}(t) (1 - \\lambda \\Delta t) \\mu \\Delta t \\\\\n",
    "&+& P_n(t) \\lambda \\Delta t \\mu \\Delta t + o(\\Delta t)\n",
    "\\end{array}\n",
    "$$\n",
    "忽略全部 $o(\\Delta t)$ 项，有\n",
    "$$\n",
    "\\frac{P_n(t + \\Delta t) - P_n(t)}{\\Delta t} = \\lambda P_{n - 1}(t) \n",
    "+ \\mu P_{n + 1}(t) - (\\lambda + \\mu)P_n(t) + \\frac{o(\\Delta t)}{\\Delta t}.\n",
    "$$\n",
    "令 $\\Delta t \\to 0$，有\n",
    "$$\n",
    "\\frac{d}{d t} P_n(t) = \\lambda P_{n - 1}(t) + \\mu P_{n + 1}(t) - (\\lambda + \\mu)P_n(t).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以上常微分方程对 $n > 0$ 成立。对 $n = 0$，因为 $n - 1$ 不存在，所以情况 (3) 无意义，且情况 (4) 对于充分小的 $\\Delta t$ 可以忽略，于是\n",
    "$$\n",
    "P_0(t + \\Delta t) = P_0(t) (1 - \\lambda \\Delta t) + P_1(t)(1 - \\lambda \\Delta t) \\mu \\Delta t,\n",
    "$$\n",
    "有\n",
    "$$\n",
    "\\frac{d}{d t}P_0(t) = -\\lambda P_0(t) + \\mu P_1(t).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当 $P_n(t)$ 为稳态解时，有 $P_n(t) = P_n$ 与 $t$ 无关，即\n",
    "$$\n",
    "\\frac{d}{d t} P_n(t) = 0, \\forall n,\n",
    "$$\n",
    "也即\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{l}\n",
    "\\lambda P_{n - 1} + \\mu P_{n + 1} - (\\lambda + \\mu) P_n = 0 \\\\\n",
    "-\\lambda P_0 + \\mu P_1 = 0.\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "而这是一个关于 $P_n$ 的差分方程，且初值满足\n",
    "$$\n",
    "P_1 = \\frac{\\lambda}{\\mu} P_0 = \\rho P_0.\n",
    "$$\n",
    "代入方程得：\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "&\\lambda P_0 + \\mu P_2 - (\\lambda + \\mu) P_1 = 0 \\\\\n",
    "\\Rightarrow & \\lambda P_0 + \\mu P_2 = (\\lambda + \\mu) \\rho P_0 \\\\\n",
    "\\Rightarrow & P_2 = \\frac{1}{\\mu}\\left(\\frac{(\\lambda + \\mu)\\lambda}{\\mu} - \\lambda\\right) P_0 = \\frac{\\lambda^2}{\\mu^2} P_0 = \\rho^2 P_0.\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一般地，有\n",
    "$$\n",
    "P_n = \\rho^n P_0.\n",
    "$$\n",
    "对 $\\rho = \\frac{\\lambda}{\\mu} < 1$，由 $P_n$ 是一个概率分布，有\n",
    "$$\n",
    "\\sum_{n = 0}^\\infty P_n = 1 \\Rightarrow \\sum_{n = 0}^\\infty P_n \n",
    "= \\sum_{n = 0}^\\infty (\\rho^n P_0) = \\frac{1}{1 - \\rho}P_0 = 1,\n",
    "$$\n",
    "也即\n",
    "$$\n",
    "P_0 = 1 - \\rho.\n",
    "$$\n",
    "因此系统的稳态分布为\n",
    "$$\n",
    "P_n = (1 - \\rho)\\rho^n, \\rho < 1.\n",
    "$$\n",
    "对于这样一个显式给出的分布，我们容易求出全部关心的（期望）指标："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**系统中平均顾客数$L_s$**\n",
    "$$\n",
    "L_s = \\sum_{n = 0}^\\infty n P_n \n",
    "= \\sum_{n = 0}^\\infty n (1 - \\rho) \\rho^n = \\frac{\\lambda}{\\mu - \\lambda}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**处于等待状态的平均顾客数$L_q$**\n",
    "\n",
    "对于 M/M/1 系统，一个顾客在接受服务的时候，其他顾客都在等。当只有一位顾客时，没有人在等，所以\n",
    "$$\n",
    "\\begin{array}{rcl}\n",
    "L_q &=& \\sum_{n = 1}^\\infty (n - 1)P_n\\\\\n",
    "&=& \\sum_{n = 1}^\\infty (n - 1)(1 - \\rho) \\rho^n \\\\\n",
    "&=& \\sum_{n = 1}^\\infty n(1 - \\rho)\\rho^n - \\sum_{n = 1}^\\infty (1 - \\rho)\\rho^n \\\\\n",
    "&=& L_s - \\rho \\\\\n",
    "&=& \\frac{\\rho \\lambda}{\\mu -\\lambda}.\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**顾客平均逗留时间$W_s$**\n",
    "\n",
    "这完全是另一个故事了，它服从\n",
    "$$\n",
    "f(w) = (\\mu - \\lambda)e^{-(\\mu - \\lambda)w}\n",
    "$$\n",
    "于是\n",
    "$$\n",
    "W_s = \\int_0^{\\infty} w f(w) d w = \\frac{1}{\\mu - \\lambda}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**平均等待时间**\n",
    "$$\n",
    "W_q = W_s - \\frac{1}{\\mu} = \\frac{\\rho}{\\mu - \\lambda}.\n",
    "$$\n",
    "这里 $\\frac{1}{\\mu}$ 的意义是平均服务时间。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "到此为止，我们只是做了一个数学的游戏，直接解析地求得了全部需要的结果。在实际工作中，我们几乎不可能做到这样完整。但是计算机模拟确可以轻松地估计上述我们想要的结果。于是，以上的解析结果，不但给我们提供了一个构思计算机模拟过程的直观引导，同时，也提供了一个用于测试模拟结果是否正确的样板。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 医院预检处的模拟程序\n",
    "\n",
    "我们来考察一个实际模拟。基本原则只有一个：简单粗暴！\n",
    "\n",
    "**参数**\n",
    "+ 每天营业时间：10 小时；\n",
    "+ 顾客到达率：65 人/小时；\n",
    "+ 服务速度：60 人/小时；\n",
    "+ 顾客总人数：营业时间 * 到达率 * 2（过宽估计）\n",
    "+ 服务模型：M/M/1\n",
    "\n",
    "注意这里 $\\rho = \\frac{\\lambda}{\\mu} = \\frac{60}{65} < 1$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "# 这句话是非标准的python，用于ipthon或jupyter这样的系统中，表示绘图即刻自动展开。\n",
    "%matplotlib inline\n",
    "\n",
    "# 这里把全部Warning过滤掉了. \n",
    "# 参见https://docs.python.org/2/library/warnings.html\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "#np.random.seed(250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_time = 10\n",
    "lam = 65\n",
    "mu = 60\n",
    "no_total_guests = total_time * lam * 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**抽样**\n",
    "\n",
    "顾客，由 $\\lambda = 65$，以及服务模型，服从指数分布，即\n",
    "$$\n",
    "f_T(t) = \\lambda e^{-\\lambda t}, t > 0,\n",
    "$$\n",
    "因此可以直接抽样生成全部顾客到达时间间隔。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arriving_gap = np.random.exponential(1 / lam, no_total_guests)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "把它们加一下，就是顾客陆续到达的时间了。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arriving_time = np.cumsum(arriving_gap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于受营业时间限制，所以就把过了营业时间的顾客都排除。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arriving_time = arriving_time[arriving_time < total_time]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "再更新一下顾客总数。注意，严格来说，这个可能仍然不是真正的当天服务顾客总数。因为要看当营业时间结束时，是否允许队列中的顾客继续等待。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_total_guests = len(arriving_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在生成每个顾客的服务时间，仍然服从指数分布："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "service_time = np.random.exponential(1 / mu, no_total_guests)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**统计**\n",
    "到目前为止，每个顾客的到达时间和服务时间，都是并行独立抽取的。事实上，整个服务系统一天的服务信息都已经确定，接下去只是统计的工作。但是这个统计是次序依赖的：\n",
    "\n",
    "+ 第一位顾客等待时间必然是 $0$；\n",
    "+ 第一位顾客离开时间必然为到达时间 + 服务时间；\n",
    "+ 假设我们已经知道第 $i - 1$ 位顾客的全部信息，那么对第 $i$ 位顾客：\n",
    "    + 等待时间 = max (第 $i - 1$ 位顾客离开时间 - 第 $i$ 位顾客到达时间, 0);\n",
    "    + 离开时间 = 到达时间 + 等待时间 + 服务时间；"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "waitting_time = np.zeros(no_total_guests)\n",
    "leaving_time = np.zeros(no_total_guests)\n",
    "for i in range(1, no_total_guests):\n",
    "    waitting_time[i] = np.max([leaving_time[i - 1] - arriving_time[i], 0])\n",
    "    leaving_time[i] = arriving_time[i] + waitting_time[i] + service_time[i]\n",
    "leaving_time[0] = arriving_time[0] + service_time[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "问一个简单的问题，在这样的系统下，我们是否应该早一点去排队？这个问题的回答似乎可以是早晚无所谓，只要抢不到头几个，接下去似乎稳定了都一样。那么模拟的结果是如何呢？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(waitting_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以看到，等待时间是分阶段递增的，也就是过一段时间，平均等待时间会上一个平台。所以结论也是显然的：越早去排队越好！"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意到我们的计算机模拟，事实上是独立于理论分析的。所以进一步可以通过统计来验证理论分析的正确性，或者说，互相验证。对于更加复杂的排队规则，计算机模拟几乎是不敏感的，然而即便对于最简单的模型，模拟的结论并不是非常显然的。因此每次对规则作出调整，适当规模的计算机模拟是必须的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后我们模拟一下队伍长度的动态变化，首先当第一个顾客到达时，队伍长度记录第一个变化时间点为第一个顾客到达时间，队伍初始长度为 $1$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "queue_length = np.zeros((2, 2 * no_total_guests))\n",
    "arriving_no = 1\n",
    "leaving_no = 0\n",
    "k = 1\n",
    "queue_length[0, 0] = 1\n",
    "queue_length[1, 0] = arriving_time[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去类似归并排序算法，在顾客到达时间序列和顾客离开时间序列中取更早的一个作为下一个队伍长度变化的时间点，如果是到达，则队伍长度加一，如果是离开，则减一。注意，考虑实际意义，到达序列一定先于离开序列结束。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "while (arriving_no < no_total_guests):\n",
    "    if (arriving_time[arriving_no] < leaving_time[leaving_no]):\n",
    "        queue_length[0, k] = queue_length[0, k - 1] + 1\n",
    "        queue_length[1, k] = arriving_time[arriving_no] \n",
    "        arriving_no += 1\n",
    "    else:\n",
    "        queue_length[0, k] = queue_length[0, k - 1] - 1\n",
    "        queue_length[1, k] = leaving_time[leaving_no]\n",
    "        leaving_no += 1\n",
    "    k += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，已经不会有新顾客到达，只需依次记录剩余顾客离开时间，队伍长度依次减一。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "while (leaving_no < no_total_guests):\n",
    "    queue_length[0, k] = queue_length[0, k - 1] - 1\n",
    "    queue_length[1, k] = leaving_time[leaving_no]\n",
    "    leaving_no += 1\n",
    "    k += 1    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "看一下实际结果："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(queue_length[1, :], queue_length[0, :])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
